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1. INTOCDOCTIDN 

Accurate determination of the thermal stresses induced in hot section 
components remains one of the most difficult problems facing engine 
design/analysts. There currently exists no rational analytical nor 
numerical techniques which can effectively deal with this problem. 
Analysts involved in hot fluid dynamics using the finite difference method 
have little interaction with those engaged in thermal stress analysis where 
the finite element method is dominant. However, the temperature 
distribution in many structural components is strongly influenced by the 
external hot gas flow, the internal cooling system of the component, and 
the structural deformation. As a result, the only effective way to deal 
with this problem is to develop an integrated solid mechanics, fluid 
mechanics, and heat transfer approach. 

In the present work, the boundary element method (BEM) is chosen as 
the basic analysis tool principally because the definition of temperature, 
flux, displacement and traction are very precise on a boundary-based 
discretization scheme. One fundamental difficulty is, of course, that a 
BEM formulation requires a considerable amount of analytical work, which is 
not needed in the other numerical methods. 

This report details progress nade, during the period November 1987 - 
November 1988 in a multi-year program commencing in March 1986, toward the 
development of a boundary element formulation for the study of hot fluid- 
structure interaction in Earth-to-Orbit engine hot section components. The 
primary thrust of the program to date has been directed quite naturally 
toward the examination of fluid flow, since boundary element methods for 
fluids are at a much less developed state. 

During the first year, work focused on the completion of a 
comprehensive literature review of integral methods in fluids, the 
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development of integral formulations for both the solid and fluid, and some 
preliminary infrastructural enhancements to a boundary element code to 
permit incorporation of the fluid-structure problem. In the second year, 
emphasis shifted to the implementation and validation phases. Boundary 
element formulations were implemented in two-dimensions for both the solid 
and the fluid. The solid was modeled as an uncoupled thermoelastic medium 
under plane strain conditions, while several formulations were investigated 
for the fluid. For example, both vorticity and primative variable 
approaches were implemented for viscous, incompressible flow, and a 
compressible version was developed. All of the above boundary element 
implementations were incorporated in a general purpose two-dimensional 
code. Thus, problems involving intricate geometry, multiple generic 
modeling regions, and arbitrary boundary conditions are all supported. 
Further details can be found in Dargush et al (1986, 1987). 

In the early portion of this past year, a number of significant 
advances were made. First, two-dimensional integration schemes were 
enhanced to obtain more accurate coefficients with somewhat less computing 
effort. This improvement was found to be particularly beneficial for 
incompressible flow, where the precise determination of the coefficients is 
imperative. Secondly, both full and modified Newton- Raphson algorithms 
were developed. This greatly inproved the convergence characteristics of 
the set of nonlinear equations governing viscous flow. Additionally, a 
region-by-region reference velocity was introduced into the formulation to 
shift the highly nonlinear portion away from the free stream and toward 
obstacles and walls, where a more refined model is appropriate. 

The combination of these advances permits the solution of a wide 
variety of thermoviscous flow problems in the low to m6derate Reynolds 
number range. Several examples are included in this report. However, at 
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higher Reynolds numbers, there is a need to get more of the physics of the 
problem into the boundary element fundamental solution. Consequently, the 
development of new convective fundamental solutions and integral 
formulations has been the primary focus of our most recent efforts. 

In the next section, a brief review of the applicable boundary element 
literature is presented. This is followed by the development of integral 
formulations for the solid in Section 3 and for the fluid in Section 4. 
Several detailed numerical examples are presented at the end of each of 
those two sections. In the fluids portion, development of the new 
convective formulations is emphasized. The remaining sections then 
summarize the progress achieved to date, and outline the work plan for the 
next year. Tables and figures appear at the end of the corresponding 
section, while references are provided in Appendix a. 

2. LITERATURE REVIEW 

Virtually nothing has appeared in the literature on the analysis of 
coupled thermoviscous fluid/ structure problems via the boundary element 
method, although some work has been done on the fluid and solid separately. 
In general, the solid portion of the problem has been addressed to a much 
greater degree. For example, a boundary-only steady-state thermoelastic 
formulation was initially presented by Cruse et al (1977) and Rizzo and 
Shippy (1977). Recently, the present authors developed and implemented the 
quasistatic counterpart (Dargush, 1987; Dargush and Banerjee, 1988a,b), 
which is presented in detail in Section 3. Others, notably Sharp and 
Crouch (1986) and Chaudouet (1987), introduce volume integrals, to 
represent the equivalent thermal body forces. A similar domain based 
approach was taken earlier by Banerjee and Butterfield (1981) in the 
context of the analogous geomechanical problem. 

An extensive review of the applications of integral formulations to 
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viscous flow problems was included in the previous annual report (Dargush. 
et al, 1987), and will not be repeated here. Interestingly, only a few 
groups of researchers are actively pursuing the further development of 
boundary elements for the analysis of viscous fluids. The work reported in 
Piva and Morino (1987) and Piva et al (1987) focuses heavily on the 
development of fundamental solutions and integral formulations with little 
emphasis on implementation. On the other hand, Tosaka and Kakuda (1986, 
1987), Tosaka and Onishi (1986) have implemented single region boundary 
element formulations using approximate incompressible fundamental 
solutions. This latter group has developed sophisticated non-linear 
solution algorithms, and consequently, are able to demonstrate moderately 
high Reynolds number solutions. Meanwhile, as will be seen in Section 4, 
the present work represents a significant advancement in the state-of-the- 
art from both a formulation and implementation standpoint. 

3. INTEGRAL FORMULATION FOR SCUDS 

3.1 Introduction 

In the current section, a surface only time domain boundary element 
method will be described for a thermoelastic body under quasistatic 
loading. Thus, transient heat conduction is included, but inertial effects 
are ignored. Formulations have been developed for three-dimensional, two- 
dimensional and axisyrrmetric problems (Dargush, 1987; Dargush and Banerjee, 
I988a,b), however, only the 2D plane strain case is detailed below. 
Separate subsections present the governing differential equations, the 
integral equations, and an overview of the numerical implementation. 

3.2 Governing Equations 

With the solid assumed to be a linear thermoelastic medium, the 
governing differential equations for transient thermoelasticity can be 
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written: 


U+4) 


9 2 u, aV 

J_ + ^ — 

3xj0Xj 


(3X+2p) a — “ 0 
dx i 


(3.1a) 


ae , a 2 © 

pc e K = k 


where 


u^ displacement vector 

© temperature 

t time 


Lagrangian coordinate 
k thermal conductivity 

p mass density 

c c specific heat at constant deformation 

X, |i Lame's constants 

o coefficient of thermal expansion 


(3.1b) 


Standard indicial notation has been employed with summations indicated 
by repeated indices. For two-dimensional problems considered herein, the 
Latin indices i and j vary from one to two. 

Note that (3.1b) is the energy equation and that (3.1a) represents the 
momentum balance in terms of displacements and temperature. The theory 
portrayed by the above set of equations, formally labeled uncoupled 
quasistatic thermoelasticity, can be derived from thermodynamic principles. 
(See Boley and Weiner (i960) for details.) 


3.3 Integral Representations 

Utilizing equation (3.1) for the solid along with a generalized form 
of the reciprocal theorem, permits one to develop the following boundary 
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integral equation: 


c„ ($)u 0 ($,t) = / tG„ *t 0 (X,t) - F. *u.(X,t)]dS(X) . (3.2) 

P<» P s pa p Pa p 

where 

a, p indices varying from l to 3 

s surface of solid 

u a> t a generalized displacement and traction 

T 

u = [u. u. 0] 
a i z 

fc a = lt l fc 2 «] T 

0,q temperature, heat flux 

G ap’ F ap generalized displacement and traction kernels (Dargush, 

1987,1988a) 

c constants determined by the relative smoothness of s at ? 


and, for example, 
t 

G ap* fc a = J G ap (x ' tj ? ’ T) fc a (x ' T> dT 
0 

denotes a Riemann convolution integral. 

In principle, at each instant of time progressing from time zero, this 
equation can be written at every point on the boundary. The collection of 
the resulting equations could then be solved simultaneously, producing 
exact values for all the unknown boundary quantities. In reality, of 
course, discretization is needed to limit this process to a finite number 
of equations and unknowns. Techniques useful for the discretization of 
(3.2) are the subject of the following sectioa 
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3.4 Numerical Implementation 

3.4.1 Introduction 

Hie boundary integral equation (3.2), developed in the last section, 
is an exact statement. No approximations have been introduced other than 
those used to formulate the boundary value problem. However, in order to 
apply (3.2) for the solution of practical engineering problems, 
approximations are required in both time and space. In this section, an 
overview of a general-purpose, state-of-the-art numerical implementation is 
presented. Many of the features and techniques to be discussed, in this 
section, were developed previously for elastostatics (e.g., Banerjee et al, 
1985,1988), and elastodynamics (e.g., Banerjee et al, 1986; Ahmad and 
Banerjee, 1988), but are here adapted for thermoelastic analysis. 

3.4.2 Temporal Discretization 

Consider, first, the time integrals represented in (3.2) as 
convolutions. Clearly, without any loss of precision, the time interval 
from zero to t can be divided into N equal increments of duration At. 

By assuming that the primary field variables, t^ and up, are constant 
within each At time increment, these quantities can be brought outside of 
the time integral. That is, 

N nAt 

• • 

Gp a *t p (X,t) » J tp(X) | Gp a (X-5.t-T)dr (3.3a) 

n=l (n-l)At 

N nAt 

F pa*V X,t) = 5 Up(X) J Fp a (X-£.t-T)dT , (3.3b) 

n=l (n-l)At 

where the superscript on the generalized tractions and displacements, 
obviously, represents the time increment number. Notice, also, that, 
within an increment, these primary field variables are now functions of 
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position only. Next, since the integrands remaining in (3.3) are known in 
explicit form from the fundamental solutions, the required temporal 
integration can be performed analytically, and written as 


N+l-n 

nAt 


G P a (X-?) 

* j Gp 0 (X-5,t-T)dr 

(3.4a) 

(n-l)At 


N+l-n 

nAt 

• 


V «-«) 

= J F pa (X-5.t-t)dt . 

(3.4b) 


(n-l)At 


These kernel functions, GjJ a (x-$) and Fj} a (x-?), are detailed in Appendix B. 
Combining (3.3) and (3.4) with (3.2) produces 

N N+l-n N+l-n 

C p O (0uJ(5) =2 J [ Gp 0 (X-5)tJ(X) - F po <X-5)uj}(X) j dS(X) , 

n=l s (3 .5) 

which is the boundary integral statement after the application of the 
temporal discretization. 


3.4.3 Spatial Dis cretization 

With the use of generalized primary variables and the incorporation of 
a piecewise constant time stepping algorithm, the boundary integral 
equation (3.5) begins to show a strong resemblance to that of 
elastostatics, particularly for the initial time step (i.e., N=l). In this 
subsection, those similarities will be exploited to develop the spatial 
discretization for the coupled quasistatic problem with two-dimensional 
geometry. This approximate spatial representation will, subsequently, 
permit numerical evaluation of the surface integrals appearing in (3.5). 
The techniques described here, actually, originated in the finite element 
literature, but were later applied to boundary elements by Lachat and 
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Watson (1976). 

The process begins by subdividing the entire surface of the body into 
individual elements of relatively simple shape. The geometry of each 
element is, then, completely defined by the coordinates of the nodal points 
and associated interpolation functions. That is, 


X( c> = x i( o * VC)x iw 


(3.6) 


with 


Nw 

X 


iw 


intrinsic coordinates 
shape functions 
nodal coordinates 


and where w is an integer varying from one to W, the number of geometric 
nodes in the element. Next, the same type of representation is used, 
within the element, to describe the primary variables. Thus, 


u S<c> - V c )uJ 


aw 


fc Su> - vc)t£ 


CUt> 


(3.7a) 


(3.7b) 


in which u£ u and t£ w are the nodal values of the generalized displacement 
and tractions, respectively, for time step n. Also, in (3.7), the integer 

w varies from one to fl, the total number of functional nodes in the 
element. Prom the above, note that the same number of nodes, and 
consequently shape functions, are not necessarily used to describe both the 
geometric and functional variations. Specifically, in the present work, 
the geometry is exclusively defined by quadratic shape functions. In two- 
dimensions, this requires the use of three-noded line elements. On the 
other hand, the variation of the primary quantities can be described, 
within an element, by either quadratic or linear shape functions. (The 
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introduction of linear variations proves computationally advantageous in 
some instances.) 

Once this spatial discretization has been accomplished and the body 
has been subdivided into N elements, the boundary integral equation can be 
rewritten as 

N M N+l-n 

C Po ( 5 )u P ( ^ = l t 5 J [ Gft a (X(C)-5)N (1) (C)tft u 
n=l m=l Sjn 


N+l-n 

- ^*p a (X( ? )- 5 ) N (d ( ?)Up w ]dS(X(c)) } 

where 


s 'I I ** ■ 


(3.8) 


In the above equation, t^ w and up^ are nodal quantities which can be 
brought outside the surface integrals. Thus, 


N M N+l-n 

c p a (5) U p(5) =5(5 t J G pa (X(c)-«)N tt (c)dS(X(c)) 

n=l m=l Sm 

N+l-n 

" U JL | Fp a (X( c)-«)N u (c)dS(X(c)) }. (3.9) 

% 

The positioning of the nodal primary variables outside the integrals is, of 
course, a key step, since now the integrands contain only known functions. 
However, before discussing the techniques used to numerically evaluate 
these integrals, a brief discussion of the singularities present in the 
kernels Gj* a and pj^ i s i n order. 

The fundamental solutions to the uncoupled quasistatic problem contain 
singularities when the load point and field point coincide, that is, when 
r=0. The same is true of Gjj a and Fp a , since these kernels are derived 
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directly from the fundamental solutions. Series expansions of terms 
present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A nurrber of observations concerning the results of these expansions 
should be mentioned. First, as would be expected, has a stronger level 
of singularity than does the corresponding G^, since an additional 
derivative is involved in obtaining F*p from G*p. Second, the coupling 
terms do not have as a high degree of singularity as do the corresponding 
non-coupling terms. Third, all of the kernel functions for the first time 
step could actually be rewritten as a sum of steady-state and transient 
components. That is. 






+ 


tr p i 
ap • 


Then, the singularity is completely contained in the steady-state portion. 
Furthermore, the singularity in G^j and F^j is precisely equal to that for 
elastostatics, while the G^ and F^ singularities are identical to those 
for potential flow. (For two-dimensions, the subscript e equals three.) 
This observation is critical in the numerical integration of the kernel 
to be discussed in the next subsection. However, from a physical 
standpoint, this means simply that, at any time t, the nearer one moves 
toward the load point, the closer the quasistatic response field 
corresponds with a steady-state field. Eventually, when the sampling and 
load points coincide, the quasistatic and steady-state responses are 
indistinguishable. As a final item, after careful examination of Appendix 
B, it is evident that the steady-state components in the kernels G^ and 

F£ p , with n>l, vanish. In that case, all that remains is a transient 
portion that contains no singularities. Thus, all singularities reside in 
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the ss G 0 p and ss F Q p components of G*p and F*p, respectively. 

3.4.4 flumeiigal integration 

Having clarified the potential singularities present in the coupled 
kernels, it is now possible to consider the evaluation of the integrals in 
equation (3.9). That is, for any element m, the integrals 

J S Gpa 1 " n ^(c)-ON (0 (C )dS(X(C)) (3.10a) 

m F 

*S FK 1 " n (X(c)-?)N M (OdS(X(c)) (3.10b) 

m * 

will be examined. To assist in this endeavor, the following three distinct 
categories can be identified: 

(1) The point $ does not lie on the element m 

(2) The point £ lies on the element m, but only non-singular or 
weakly singular integrals are involved 

(3) The point $ lies on the element m, and the integral is strongly 
singular. 

In practical problems involving many elements, it is evident that most 
of the integration occurring in equation (3.9) will be of the Category (1) 
variety. In this case, the integrand is always non-singular, and standard 
Gaussian quadrature formulas can be employed. Sophisticated error control 
routines are needed, however, to minimize the computational effort for a 
certain level of accuracy. This non-singular integration is the most 
expensive part of a boundary element analysis, and, consequently, must be 
optimized to achieve an efficient solution. In the present implementation, 
error estimates, based upon the work of Stroud and Secrest (1966), are 
employed to automatically select the proper order of the quadrature rule. 
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Additionally, to improve accuracy in a cost-effective manner, a graded 
subdivision of the element is incorporated, especially when ? is nearby. 
For two-dimensional problems, the integration order varies from two to 
twelve, within each of up to four element subdivisions. 

Turning next to Category (2), one finds that again Gaussian quadrature 
is applicable, however, a somewhat modified scheme must be utilized to 
evaluate the weakly singular integrals. This is accomplished in two- 
dimensional elements via suitable subsegmentation along the length of the 
element so that the product of shape function, Jacobian and kernel remains 
well behaved. 

Unfortunately, the remaining strongly singular integrals of Category 
(3) exist only in the Cauchy principal value sense and cannot, in general, 
be evaluated numerically, with sufficient precision. It should be noted 
that this apparent stumbling block is limited to the strongly singular 

portions, ss F^j and SSf 00 , of the F*p kernel. The remainder of F*p, 
including tr F^ and tr F 0e » can be computed using the procedures outlined 
for Category (2). However, as will be discussed in the next subsection, 
even the Category (3) ss F^j and ss p 00 kernels can be accurately determined 
by employing an indirect 'rigid body' method originally developed by Cruse 
(1974). 


3.4.5 As.sentoly 

The complete discretization of the boundary integral equation, in both 
time and space, has been described, along with the techniques required for 
numerical integration of the kernels. Now, a system of algebraic equations 
can be developed to permit the approximate solution of the original 
quasistatic problem. This is accomplished by systematically writing (3.9) 
at each global boundary node. The ensuing nodal collocation process, then. 
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produces a global set of equations of the form 
N 

J ( tG N+1 " n Ht n } - [F N+1 ~ n l tu n J ) = (0} , (3.11) 

n=l 

where 

[G N+1_n ] unassembled iratrix of size (d+l)P x (d+l)Q, with 
coefficients determined from (3.10a) 

[pto+l-nj assembled matrix of size (d+l)P x (d+l)P, with coefficients 
determined from (3.10b) and c p<x included in the diagonal 
blocks 

ft n } global generalized nodal traction vector with (d+l)Q 

components 

(u n ) global generalized nodal displacement vector with (d+l)P 
components 

{0} null vector with (d+l)P components 

P total number of global functional nodes 

M 

Q - K 

m=l 

A m number of functional nodes in element m 

d dimensionality of the problem. 

In the above, recall that the terms generalized displacement and traction 
refer to the inclusion of the temperature and flux, respectively, as the 
(d+l) component at any point. 

Consider, now, the first time step. Thus, for N=i, equation (3.11) 
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becomes 


IG 1 ] tt 1 } - IF 1 ! tu 1 ) = CO) . (3.12) 

However, at this point, the diagonal block of IF 1 ! has not been completely 
determined due to the strongly singular nature of ss F^j and ss F ee . 
Following Cruse (1974) and, later, Banerjee et al (1986) in elastodynamics, 
these diagonal contributions can be calculated indirectly by imposing a 
uniform 'rigid body' generalized displacement field on the same body, but 
under steady-state conditions. Then, obviously, the generalized tractions 
must be zero, and 

I SS F] (l) = (0) , (3.13) 

where Cl) is a vector having all (d+l)P components equal to one. Using 

(3.13) . the desired diagonal blocks, ss F i; . and ss F ee , can be obtained from 
the summation of the off-diagonal terms of C SS F). The remaining transient 
portion of the diagonal block is non-singular, and hence can be evaluated 
to ary desired precision. With that step completed, (3.12) is rewritten as 

CG 1 ]{t 1 ) - IF 1 ) Cu 1 ) = CO) . (3.14) 

In a well-posed problem, at time At, the set of global generalized 
nodal displacements and tractions will contain exactly (d+DP unknown 
components. Then, as the final stage in the assembly process, equation 

(3.14) can be rearranged to form 

IA 1 ! (x 1 ) = IB 1 ) Cy 1 ) , (3.15) 

in which 

Cx 1 ) unknown components of Cu 1 ) and Ct 1 ) 

Cy 1 } known components of Cu 1 ) and Ct 1 ) 
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[A 1 ]. IB 1 ! associated coefficient matrices. 


3.4.6 Solution 

To obtain a solution of (3.15) for the unknown nodal quantities, a 
decomposition of matrix [A 1 ] is required. In general, [A 1 ] is a densely 
populated, un symmetric matrix. The out-of-core solver, utilized here, was 
developed originally for elastostatics from the L INPACK software package 
(Dongarra et al, 1979) and operates on a submatrix level. Within each 
submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [A 1 ] is stored in a 
direct-access file for reuse in subsequent time steps. Backsubstitution 
then completes the determination of Cx 1 }. Additional information on this 
solver is available in Banerjee et al (1985). 

After returning from the solver routines, the entire nodal response 
vectors, {u 1 } and (t 1 ), at time At are known. For solutions at later 
times, a simple marching algorithm is employed. Thus, from (3.11) with 
N=2, 

[G 2 ] (t 1 ) - [F 2 ] Cu 1 J + IG 1 ] (t 2 ) - IF 1 1 Cu 2 ) = {0} . (3.16) 

Assuming that the same set of nodal components are unknown as in (3.14) for 
the first time step, equation (3.16) is reformulated as 

IA 1 Hx 2 ) = [B 1 ] Cy 2 ) - [G 2 ] (t 1 ) + IF 2 ] (u 1 ) . (3.17) 

Since, at this point, the right-hand side contains only known quantities, 
(3.17) can be solved for (x 2 ). However, the decomposed form of [A 2 ] 
already exists on a direct-access file, so only the relatively inexpensive 
backsubstitution phase is required for the solution. 

The generalization of (3.17) to any time step N is simply 
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(3.18) 


N-l 

[A 1 ]^} = [B 1 ]^} - J ( [G N+1 ' n Ht n } - [F N+1 " n ] (U n ) ) 

n=i 

in which the summation represents the effect of past events. By 
systematically storing all of the matrices and nodal response vectors 
computed during the marching process, surprisingly little computing time is 
required at each new time step. In fact, for any time step beyond the 
first, the only major computational task is the integration needed to form 
[G n ] and [F n 1. Even this process is somewhat simplified, since now the 
kernels are non-singular. Also, as time marches on, the effect of events 
that occurred during the first time step diminishes. Consequently, the 
terms containing IG N ] and [F N ] will eventually become insignificant 
compared to those associated with recent events. Once that point is 
reached, further integration is unnecessary, and a significant reduction in 
the computing effort per time step can be achieved. 

It should be emphasized that the entire boundary element method 
developed, in this section, has involved surface quantities exclusively. A 
complete solution to the well-posed linear uncoupled quasistatic problem, 
with homogeneous properties, can be obtained in terms of the nodal response 
vectors, without the need for ary volume discretization. In many practical 
situations, however, additional information, such as, the temperature at 
interior locations or the stress at points on the boundary, is required. 
The next subsection discusses the calculation of these quantities. 

3.4.7 nr Qia ptitieS 

Once equation (3.18) is solved, at any time step, the complete set of 
primary nodal quantities, (u N ) and Ct N l, is known. Subsequently, the 
response at points within the body can be calculated in a straightforward 
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manner. 


For any point 5 in the interior, the generalized displacement can 
be determined from (3.9) with c pa = 6 pa . That is, 

N M 

u 2<{> - 5 1 5 1 tj. /s Gg; 1 ' n <x(t)-{)N <J (c)ds(x(o) 
n=l m=l m 

- ujj f s F^ 1 ' n (X(c)-?)N u (0dS(X(0) ]) . (3.19) 

p & m p 

Now, all the nodal variables on the right-hand side are known, and. as long 
as, 5 is not on the boundary, the kernel functions in (3.19) remain non- 
singular. However, when ? is on the boundary, the strong singularity in 
ss 

F pa prohibits accurate evaluation of the generalized displacement via 
(3.19), and an alternate approach is required. The apparent dilemma is 
easily resolved by recalling that the variation of surface quantities is 
completely defined by the elemental shape functions. Thus, for boundary 
points, the desired relationship is simply 

U 2u> = VO U 1 (3 * 20) 

where N w ( are the shape functions for the appropriate element and 
c are the intrinsic coordinates corresponding to $ within that element. 
Obviously, from (3.20), neither integration nor the explicit contribution 
of past events are needed to evaluate generalized boundary displacements. 

In many problems, additional quantities, such as heat flux and stress, 
are also important. The boundary integral equation for heat flux, can be 
written 

N M 

- 5 ( l I / s E^ 1 1 ' n (X(c)-?m„«><3S<X(i)> 
n=l m=l m 

- un w f s D JJi" n (X(c)-ON (i> (c)dS(X(c)) 1) (3.2D 
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where 


E p©i<*( 0-5> - - k H 


9Gp fl (X (c)-C) 
i 


(3.2la) 


0F g©(X(O-?) 
= - k g^T 


(3.2lb) 


This is valid for interior points, whereas, when 5 is on the boundary, the 
shape functions can again be used. In this latter case. 


V c >qE = 


dN » (C) qN _ I qN, ,v 

3 c to " " k 3 ; q i (?) ’ 


(3.22a) 


(3.22b) 


which can be solved for boundary flux. Meanwhile, interior stresses can be 
evaluated from 



N M 


- J < 1 l t(L / s E^‘ n (X(t)-{)N 0 (OdS(X(5)) 
n=i m=l m 


- / s D p fj 1_n (X( ;)-«)V c )dS(X(c)) I) 


in which 


Ej lj (X(c)-«) 


2n v 

6 ij 


8G H t 

*r* 




f aG Pi 


P 6 ijG?© 


D p ij (x(c)-5) 


2fi v 

T=TZ 



a ?j ®«i 


p8 ij F ?« • 


(3.23) 


(3.23a) 


(3.23b) 


Equation (3.23) is, of course, developed from (3.19). Since strong kernel 
singularities appear when (3.23) is written for boundary points, an 
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alternate procedure is needed to determine surface stress. This alternate 
scheme exploits the interrelationships between generalized displacement, 
traction, and stress and is the straightforward extension of the technique 
typically used in elastostatic implementations (Cruse and Van Buren, 1971). 
Specifically, the following can be obtained 

= VOtJ. <3. 24a) 
ff ijU> - <u k.l < « ,+u i.lc < * )) * " P 8 ijV 5)u <L (3.24b) 


!!i 

3c 


u i.j<5> 


3N,, 

3 C 



(3.24c) 


in which u^ (Piously the nodal temperatures, and, 

^jkl = X8 ij 5 kl + 2 ^ 5 ik 6 jl • 

Equations (3.24) form an independent set that can be solved numerically for 
v^jU) and u^ j(?) completely in terms of known nodal quantities u^ and 
t^ w , without the need for kernel integration nor convolution. Notice, 
however, that shape function derivatives appear in (3.24c), thus 
constraining the representation of stress on the surface element to 
something less than full quadratic variation. The interior stress kernel 
functions, defined by (3.23), are also detailed in Appendix B. 

3.4.8 Advanced Features 

The thermoelastic formulation has been implemented as a segment of the 
state-of-the-art, general purpose boundary element computer program, GP- 
BEST. Consequently, many additional features, beyond those detailed above, 
are available for the analysis of complex engineering problems. Perhaps, 
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the most significant of these items, is the capability to analyze 
substructured problems. This, not only extends the analysis to bodies 
composed of several different materials, but also often provides 
corrputational efficiencies. An individual substructure or generic modeling 
region (GMR) must contain a single material. During the integration 
process, each GMR remains a separate entity. The GMR's are then brought 
together at the assembly stage, where compatibility relationships are 
enforced on cannon boundaries between regions. Typically, compatibility 
ensures continuous displacement and temperature fields across an interface, 
however, recent enhancements to the code permit sliding between regions, 
spring contacts and interfacial thermal resistance to model air gaps or 
coating resistances. In the latter instances, discontinuities appear at 
the interface. In any case, the multi-GMR assembly process produces block- 
banded system matrices that are solved in an efficient manner. 

As another feature, a high degree of flexibility is provided for the 
specification of boundary conditions. In general, time-dependent values 
can be defined in either global or local coordinates. Not only can 
generalized displacements and tractions be specified, but also spring and 
convection boundary conditions area available. Another recent addition 
permits time-dependent ambient temperatures. A final item, worthy of note, 
is the availability of a comprehensive symmetry capability which includes 
provisions for both planar and cyclic symmetry. 

These advanced features greatly extend the range of applicability of 
the present formulation. In the next section, several examples are 
presented to demonstrate the validity and applicability of this boundary- 
only formulation. 
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3.5 Examples 

3.5.1 Sudden Heat ing of an Aluminum Block 

As a first example, transient heating of an aluminum block is examined 
under plane strain conditions. The block, shown in Figure 3.1, initially 
rests in thermodynamic equilibrium at zero temperature. Then, suddenly, 
the face at Y = 1.0 in. is elevated to 100°F, while the remaining three 
faces are insulated and restrained against normal displacements. Thus, 
only axial deformation in the Y-direction is permitted. Naturally, as the 
diffusive process progresses, temperature builds along with the lateral 

stresses <j xx and <y 22 . To complete the specification of the problem, the 
following standard set of material properties are used to characterize the 
aluminum: 

E = 10x10^ psi , v = 0.33 , 

a = 13 x10 _6 /°F , 

k = 25 in.-lb./sec. in.°F , pc g = 200 in.-lb./in. 3o F . 

The two-dimensional boundary element idealization consists of the 
simple four element, eight node model included in Figure 3.1. A time step 
of 0.4 sec. is selected, corresponding to a non-dimensional time step of 
0.05. Additionally, a finite element analysis of this same problem was 
conducted using a modified thermal version of the computer code CRISP (Gunn 
and Britto, 1984). The finite element model is also a two-dimensional 
plane strain representation, however sixteen linear strain quadralaterals 
are placed along the diffusion length. In the FE run, a time step of 0.2 
sec. is employed. 

Temperatures, displacements, and stresses are compared in Table 3.1. 
Notice that the boundary element analysis, with only one element in the 
flow direction, produces a better time-temperature history than does a 
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sixteen element FE analysis with a smaller time step. Both methods exhibit 
greatest error during the initial stages of the process. This is the 
result of the imposition of a sudden temperature change. Meanwhile, the 
comparison of the overall axial displacement indicates agreement to within 
3% for the BE analysis and 5% for the FE run. A steady-state analysis via 
both methods produces the exact answer to three digit accuracy. The last 
comparison, in the table, involves lateral stresses at an integration point 
in the FE model. The boundary element results are quite good throughout 
the range, however, the FE stresses exhibit considerable error, 
particularly during the initial four seconds. Actually, these finite 
element stress variations are not unexpected in light of the errors present 
in the temperature and displacement response. Recall that in the standard 
finite element process, stresses are computed on the basis of numerical 
differentiation of the displacements, whereas in boundary elements, the 
stresses at interior points are obtained directly from a discretized 
version of an exact integral equation. Consequently, the BE interior 
stress solution more nearly coincides with the actual response. 

3.5.2 Circular Disc 

Next, transient thermal stresses in a circular disc are investigated. 
The disc of radius 'a' initially rests at zero uniform temperature. The 
top and bottom surfaces are thermally insulated, and all boundaries are 
completely free of mechanical constraint. Then, suddenly, at time zero, 
the temperature of the entire outer edge (i.e., r=a) is elevated to unity 
and, subsequently, maintained at that level. 

The boundary element model of the disc with unit radius is shown in 
Figure 3.2. Only four quadratic elements are employed, along with quarter 
symmetry. Ten interior points are also included strictly to monitor 
response. In addition, the following non-dimensional ized material 
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properties are arbitrarily selected for the plane stress analysis: 


E = 1.333 pc = i.o 
V = 0.333 k = 1.0 
a = 0.75 


Results obtained under quasistatic conditions for a time step of 0.005 are 
compared, in Figures 3.3, 3.4 and 3.5, to the analytical solution presented 
in Timoshenko and Goodier (1970). Notice that temperatures, as well as 
radial and tangential stresses are accurately determined via the boundary 
element analysis. In particular from Figure 3.5, even the tangential 
stress on the outer edge is faithfully reproduced. 

3.5.3 Turbine Blade 

For the final application, the plane strain response of an internally 
cooled turbine blade is examined under startup thermal transients. The 
boundary element model of the blade is illustrated in Figure 3.6. In this 
problem, the two GMR approach is chosen solely to enhance computational 
efficiency. This is accomplished by reducing the aspect ratio of 
individual GMR's and by creating a block banded system matrix. The leading 
(lefthand) GMR consists of 26 quadratic elements, while 24 elements are 
used to model the trailing (righthand) region. 

The blade is manufactured of stainless steel with the following 
the rmomechanical properties: 

E = 29.0 X l0 6 psi pc g = 368 in.-lb./in. 3o F 

v = 0.30 k = 1.65 in.-lb./sec.in.°F 

a = 9.6 X lO- 6 in./in.°F 

During operation a hot gas flews outside the blade, while a relatively cool 
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gas passes through the internal holes. Hie gas temperature transients are 
plotted in Figure 3.7 for a typical startup. Convection film coefficients 
are specified as follows: 

Outer surface at leading edge h = 50 in.-lb./sec. in. 2o F 
Remainder of outer surface h = 20 in.-lb./sec. in. 2o F 
Inner cooling hole surfaces h = 10 in. -lb. /sec. in. 2o F 

A time step of 0.2 sec. is employed for the boundary element analysis. 

The response at two points, A, on the leading edge and, B, at midspan 
are displayed in Figures 3.8 and 3.9. Notice that temperatures and 
stresses are consistently higher on the leading edge, reaching peak values 
of approximately 1500°F and -60 ksi, respectively. Also, as is evident 
from Figure 3.9, significant stress reversals occur during this startup. 
As a next step, these numerical results could be used as input for a 
fatigue analysis to assess the durability of the design. In that regard, 
it should be emphasized that the stresses presented for points A and B are 
surface stresses, calculated by satisfying the constitutive laws, strain- 
displacement and equilibrium directly at the boundary point. This can be 
expected to produce much more accurate results than the standard practice 
utilized in finite element approaches of extrapolating interior Gauss point 
stress values to the boundary. 


25 


TABLE 3.1 


00 



E-* 

S3 


cs 

r- 


00 

Os 

r- 

m 


CO 

• 

• 

• 

• 

• 

• 

• 

• 



v-> 

On 

^4 

CO 


wo 

VO 

r- 

•H 

1 

1 

H 

1-4 

i-4 

H 

i-4 

1-4 

CO 

w cs 

o 



1 

1 

1 

1 

1 

1 


CO co 
co wo 


a» • 

Ov 

r- 

co 

CS 

00 

o 

On 

r- 

U o w 
8 ,, S 

• 

• 

• 

• 

• 

• 

• 

• 

cn 

r- 

o 

cs 



vo 

vo 

i 

I 

iH 

iH 

Y 

^4 

r4 

^4 

M .p 



t 

1 


1 

1 

1 

!"« 

vo 

i-4 

co 

H 


WO 

CO 

O 

3 S 

• 

vo 

• 

Ov 

• 

r4 

• 

<n 

• 

• 

»o 

• 

VO 

• 

r- 

X 

l 

l 

^4 

i-4 

r4 

i-4 

1-4 

r4 

w 



1 

1 

1 

1 

1 

1 


cs 


I 


— £h 


c 

fil 

o 

O 

o 

o 

o 

o 

o 

o 

o 

•H 

CQ 

cs 

cs 

i-4 

4“ 

o 

CO 

CO 

iH 

t- 

3- 


On 

CO 

VO 

00 

o 

1-4 

N 

co 

co 



r4 

1-4 

iH 

cs 

CS 

<S 

cs 

M 


c o 

O II 

© 

© 

© 

O 

o 

© 

© 

© 

© 

as w 

VO 

vo 


vo 

CO 

r- 

00 

r- 

to 

rH >4 E 

00 

C4 

•o 

r- 

On 

© 

1-4 

cs 

CO 



iH 

H 

1-4 

1-4 

CS 

<s 

cs 

cs 


4J 

(Q 


rH 

<0 

4J 

O 

© 

© 

© 

© 

© 

8 

© 

© 

© 

•H 

<0 

1-4 

On 


00 

WO 

© 

00 


3 

X 

On 

CS 

wo 

r- 

On 

© 

CS 

CS 

co 

W 


i-4 

»-4 

1-4 

1-4 

CS 

CS 

CS 

CS 



00 

r- 

r- 

*n 

CS 

wo 

On 

OV 

00 

to 

© 

r- 

i-4 

CS 

© 

VO 

r4 

WO 


CS 

co 

WO 

vo 

C- 

l> 

00 

00 


On 


V© 

co 


o 

v*> 


O 

VO 


ON 

VO 


o 

00 


*r 

00 


4J 


© 

CO 

wo 

On 

1-4 

WO 

wo 

WO 

U 

• 

• 

• 

• 

• 

• 

• 

• 

■ 

CCS 

'•r 

CS 

00 

1-4 

1-4 

O 

vo 

1-4 

wo 

X 

w 


CS 

CO 

wo 

vo 

f"* 


00 

00 


§ v 

00 

vo 


CS 

© 

•H w 

• 

• 

• 

• 

• 

E4 w 

© 

iH 

CS 

tn 



» 'O t « 

• • t • 

^ vo r- 


26 


8.0 88.6 88.2 88.8 2400 2390 2410 -17.9 -17.7 -18.1 



RLUMINUM BLOCK 


































4 INTEGRAL FORMULATIONS FOR FLUIDS 

4.1 Introduction 

Next, attention turns to the hot fluid. In the following, a number of 
integral formulations are developed for compressible and incompressible 
thermoviscous flew and, additionally, for the simpler theory of convective 
heat transfer. Subsections present the governing equations, fundamental 
solutions, integral representations, an overview of the numerical 
implementation, a brief description of the approach for coupling the fluid 
with the solid, and, finally, a number of detailed numerical examples. 

4.2 Governing Equations 

4.2.1 Compressible Thermoviscous Flow 

The governing equations for a thermally-sensitive, compressible, 
viscous fluid can be developed from the consideration of the conservation 
laws of mass, momentum, and energy. In each case, the law is first written 
for a continuum which is, in general, moving non-uniformly with respect to 
the observer. The local (differential) form of the law is then derived. 
Although a derivation of the governing equations of fluid dynamics, similar 
to the following, can be found in a number of texts, it is a useful means 
for establishing the underlying assumptions and limitations. 

The Principle of the Conservation of Mass asserts that the time rate 
of change of mass must equal the rate of mass increase due to internal 
sources. That is, 

J pdV = J *dV , (4.1) 

V(t) V(t) 

where p is the mass density, * is the mass source rate per unit volume, and 
the operator D/Dt represents a material time derivative. Notice that in 
(4.1) the irass of interest occupies V(t), a region of space which may vary 
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with time. Applying a generalized version of Leibnitz's Rule to the left- 
hand side of (4.1) produces 

Dti pdv “I ^ + j Wfjte * (4 - 2) 

V(t) V(t) S(t) 

where S(t) is the surface enclosing V(t). and v^ and nj are the local 
velocities and outward normals on that surface, respectively. However, via 
the Divergence Theorem, the surface integral can be rewritten as 

| pv^n^dS = J fx. (pv j )<3V * <4 ’ 3) 

S(t) V(t) 3 

Therefore, from (4.1), (4.2), and (4.3) 

J [ ft + if: (pv j } " +1 ™ “ 0 * (4 * 4) 

V(t) 3 

Since this integral must vanish for all regions V(t), the integrand must be 
identically zero. Thus, 


If + afc - * - o ■ «-’> 

which is the desired local form of continuity or Conservation of Mass. 
This can also be written 


D& 

Dt 


+ 



= o , 


where 


(4.6) 


D a a 

Dt = at + v j axj 


(4.7) 
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is again the material time derivative. 

Next, consideration is given to the Conservation of Linear Momentum. 
In this case, according to Newton's Second Law, it is postulated that the 
time rate of change of momentum is equal to the resultant of the applied 
forces. Alternatively, these applied forces can be visualized as the rate 
of momentum entering the region through the surface plus the rate of 
momentum increase due to internal generation. With either interpretation, 

^ J pv £ dV = J c^njdS + J (f j+Vj 4)dV , (4.8) 

V(t) S(t) V(t) 

where <r^j is the total stress tensor and f^ is the body force vector. 
Notice that the term v^ is included in the last volume integral to 
account for the internal momentum generation due to mass sources. Applying 
the generalized Leibnitz's Rule and the Divergence Theorem to the left-hand 
integral of (4.8) yields 


Dt I pv i av - 1 f It (P V * it < ‘ ,v i v j ) 1 ® • 


V(t) 


V(t) 


(4.9) 


The Divergence Theorem can also be invoked to convert the surface integral 
in (4.8) into a volume integral. Thus, 


| . ljnj as - J sf « • 


S(t) 


3x i 
V(t) 3 


(4.10) 


Utilizing (4.9) and (4.10), Newton's Second Law becomes 


da 

at ' r 'i' dx. '’"I’D' ix ‘ 

V(t) 3 3 


I [ k (pv i } + ixT (pv i v j ) - 3^ " f i - v i* 1 " 0 • (4 - n) 


Again, since this integral must vanish for arbitrary regions, the integrand 
must be zero. That is. 
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fc ‘"V * 8x7 ,pv i v j> - ST 1 - f i - v i* ■ 0 • 

However, equation (4.12) can be rewritten as 


(4.12) 


3v. av. a 0ii 

p 3T * pv j S7 ' 5^ - f i * 


v i [ It - <pv j> - 


(4.13) 


But since the bracketed term multiplying in (4.13) equals zero from the 
continuity equation (4.5), the local form of the Conservation of Linear 
Momentum becomes 


8Vj 

W 


+ pv 


3v^ 

i5 i 



or simply 


(4.14) 


p 





(4.15) 


Note that although continuity is invoked above, a flow field that conserves 
linear momentum does not automatically conserve mass. In addition, the 
moment of momentum must also be conserved as a consequence of Newton's 
Second Law. However, satisfaction of this law only necessitates that the 
stress tensor be symmetric. 

Finally, the Conservation of Energy is examined. For energy balance, 
the time rate of change of kinetic plus internal energy must equate the 
rate of work done fcy the body forces and surface tractions, along with the 
rate of energy entering via heat transfer across the surface, the rate of 
kinetic and internal energy increase due to mass sources, and the rate of 
energy input due to heat sources. In equation form. 


- f 

Dt J 


P( 


v i v i 


+ E)dV = | f jVjdV + | w^jnjVji 


dS 


V(t) 


V(t) 


S(t) 


- j 

S(t) 
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(4.16) 


f v i v i f 

+ I i, (-“ + E)dV +1 (J> dV 

V(t) V(t) 

where E is the internal energy per unit mass, is the heat flux vector 

and 4 is the heat source rate per unit volume. By first applying the 

generalized Leibnitz's Rule to the left-hand side of (4.16), and then 

invoking the Divergence Theorem for all of the remaining surface integrals, 

equation (4.16) is transformed into 

r r n v -5 v i a 

J [ ' Bt '-p * E) - - v i £ i - ST '•ijV * 5^ - * 

V(t) 3 1 


* ( -2^ * EM ft * 8^ <pV i ) -♦ > ] - 0 • 


(4.17) 


Since this is valid for any region V(t), 

n ViVi a dq. 

p gj <— * E) - pv l £ i - 5J- + 5^ - * 

♦ <-£l ♦ E)(|E * j^Cpvj) -» ) - o . 

After further rearrangement this becomes, 


(4.18) 


DE 


3q i 


3v^ 

8x i " 

W ij 


,Vi 

2 

+ E) 

r 3p 

L at 


Dv 4 


3a. . 


0. 


(4.19) 


Now, the first bracketed expression in (4.19) vanishes via the Conservation 
of Linear Momentum, while the second bracketed expression is zero from the 
Conservation of Mass. Thus, equation (4.19) reduces to 
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(4.20) 


DE 3 ^i ”i i 

p Dt + “ °ij 3Xj " ^ " 0 


3v. 


as the expression for the Conservation of Energy. 

To recapitulate, the conservation laws for a thermoviscous fluid can 
be written collectively as 


Dt 


+ 



= 0 


Dv i _a^. 

p Dt axj E i 


= o 


Mass 


(4.2ia) 


Momentum 


(4.21b) 



3< ?i 3v i . 

"ij axj “ * 


3x 


0 . 


Energy 


(4.21c) 


Next, constitutive relationships are introduced. In particular, a 
homogeneous isotropic Newtonian fluid is assumed such that 


3v. av^ av k 

'U * 2 “ ( a^ * K? + u ij 5^ 


- 8 


ij P 


(4.22) 


where p is the thermodynamic pressure, while m and X are coefficients of 
viscosity. Fourier's law of heat conduction is also envoked, which for an 
isotropic medium becomes 



where e is the thermodynamic temperature and k is the thermal conductivity. 
Additionally, the fluid is modeled as a perfect gas. Thus, the kinetic 
equation of state is simply 


P = pRe , 


(4.24) 
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in which R is the gas constant. Finally, a relationship is needed for the 
internal energy E. From thermodynamic considerations for a perfect gas. 


h = E+ E = E + R0 (4.25a) 

P 

where h is the enthalpy. In addition, if the specific heat at constant 
pressure, c_» does not vary with temperature, then 

h = c p e . (4.25b) 

and, hence, 

E = (c p -R)e - c v e , (4.26) 

where c v is the specific heat at constant volume. Equations (4.22), 
(4.23), (4.24), and (4.26) lead to the following form of the governing 
equations for the idealized thermoviscous fluid: 


D£ 

Dt 


+ 



= 0 


p 



a 2 v 


(X+»i) 


2 - _ 


dx j* x i 


jhi + _ 

dXjSXj dx i 


f i = 


(4.27a) 


(4.27b) 


do , a 2 e 

PC,, ri*- ” k 


v Dt ax^ax^ 


where a is the viscous 


3v. 

4 = T ij a^T ' 


and the fluid stresses 


3v, 

+ p £xT“ a_( * = 0 , 


dissipation defined fcy 


(4.27c) 


(4.28) 


38 



(4.29) 


av_i 

c i 


3v b 


JVj 

T ij ■ <8^7 * + l6 lj 8^ 


Equations (4.27), along with (4.24), define a highly non-linear set of six 
equations in the six variables: velocity (v^. pressure (p), temperature 

(0), and density (p). 


4.2.2 Incompressible Thermoviscous Flow 

For incompressible flow, a number of simplifications are in order. In 
particular, the divergence of the velocity is zero, which from continuity 
requires that the density remain constant. As a result, the governing 
equations reduce to the following: 


% >\ to . 

u — : — + ~ f • 


Dt 


11 9Xj3Xj 


a 2 e 


pc v Dt " k ax.ax, 




- a - 


(4.30a) 


(4.30b) 


where 


T ij 


dv. av., 

■ 2 “ ‘SC * 


(4.31) 


3 '“I 

and the viscous dissiptation d is again defined by (4.28). It should be 
noted that now the quantity p. appearing in (4.30a), is no longer the 
thermodynamic pressure determined from (4.24), but rather the mean fluid 
pressure. 


4.2.3 Incompressible Viscous Flow 

With the assumption of isothermal conditions, the energy equation 
(4.30b) is no longer required. All that remains is the familiar Navier- 
Stokes equation 



a 2 v. 


* X j dx j 


+ ^- 

axi 



0 . 


(4.32) 
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4.2.4 fo nvec t iye. jfeat.Transfer 

On the other hand, if the flew field is known or can be approximated, 
then equation (4.30a) is superfluous. Consequently, fluid temperatures can 
be determined directly from the scalar convective-diffusion equation 


, DO _ 3 2 

v Dt k ax i 3x i 


0 . 


(4.33) 


In (4.33), the effects of viscous dissipation are included as body heat 
sources. 


4.3 Fundamental Solutions 

4.3.1 Compressible Thermoviscous Flow 

One of the primary requirements for developing a boundary element 
formulation is that the fundamental solution of the governing differential 
equations must exist. These fundamental solutions can be viewed in same 
sense as the shape functions in the finite element method. For solid 
mechanics these have been very well explored. Starting with Kelvin's 
solution (1846), investigators such as Stokes, Poisson, Boussinesq, 
Mindlin, and Nowacki have provided both static and transient solutions 
which form the basis of the boundary element formulations in solid 
mechanics. It is unfortunate that workers in fluid mechanics have not 
found much use for these fundamental solutions in the infinite space and 
therefore have not derived the corresponding fluid solutions. The 
exception is the time -dependent fundamental solution for viscous, 
incompressible Stokes flow presented in Ladyzhenskaya (1969). Since the 
boundary element formulations could not be developed without these 
solutions, a substantial amount of effort has been devoted in the present 
work to successively derive more complete solutions of the differential 
equations. In essence, each advancement brings more of the physics of the 
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problem into the fundamental solution. Below is an overview of the 
derivation for compressible, thermoviscous flow. 

As a starting point, reference values for each of the primary 
variables are introduced in an effort to produce a linearized differential 
operator. Thus, let 


v i - 0^ + Uj 

(4.34a) 

p ■ p 0 + Pa 

(4.34b) 

0 - ®o * «A 

(4.34c) 

" - P 0 * "A • 

(4.34d) 


in which U^, p Q , e Q , and p Q are constant reference values, and uj, p A » 
and p A are the perturbations. Plugging (4.34) into (4.27) yields, after 
some manipulation, 


D ,J?a 3ua 3p A D# a dUa 

M + p o 3^ - - u i * pR Dt 4 ' Pi 


(4.35a) 


D o u i a2u i 

P 0 Dt (X+p) 3 x,3xj 


d^U, 




3Uj 


D„U 


- m + m - pu . 

ijoxj dXjdXj 3x^ 3 dXj 


o u i 
P A Dt 


+ fa 


(4.35b) 


p O C v Dt K a x i 3x i 


36, 


Va 


pC v u i 3x i " P A c v Dt 


3u. 


- P + « + » 


(4.35C) 


where 


^ = L. + U . J_ 

Dt 3t i 3x t * 


(4.36) 


Now, in (4.35), the entire left-hand side involves a linear differential 
operator with constant coefficients. Notice that in the above form, the 
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operator for the energy equation involves only temperature (© A ), while the 
mass and momentum balance operators are coupled by the inclusion of both 
velocity (u^) and pressure (p A ). Terms on the right-hand side of (4.3 5) 
are, in general, non-linear, and can for the present be considered as body 
sources and forces of unknown magnitude. Then, the governing equations 
become 


’“I 7 

Dt * P 0 SIT " * 


(4.37a) 


„ 

P 0 Dt 


- (X+|i) 


a V*j ' " 


3P A _ 

+ = f , 

0Xj3Xj 3X^ i 


(4.37b) 


a 2 e, 


da. 

O A , 
p o c v Dt " K ax.axj 


a 


(4.37c) 


A fundamental solution of (4.37) is required for the boundary element 
formulation. This will be obtained subsequently, and referred to as the 
convective fundamental solution for compressible, thermoviscous flew, since 
a linearized portion of the convective derivatives are included in the 
differential operator. Interestingly, it may also be viewed as the 
fundamental solution due to stationary point forces and sources in a 
uniformly moving medium or, equivalently, as a uniformly moving point force 
and source solution in a stationary medium. The concept of moving media 
fundamental solutions is clearly developed in the excellent monograph on 
aeroacoustics by Goldstein (1976). 

Consider, first, the coupled set of equations (4.37a) and (4.37b), and 
introduce the Hemholtz decomposition of the velocity and body force, such 
that 
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(4.38a) 


Ui 



aw R 
e ijk ixT 


with 


aw. 

dXi 


i _ 


= 0 


p 


af 


aF L 


ax^ + e ijk axj 


wlth ^ ■ 0 


Ihen, (4.37b) becomes 


(4.38b) 


^I p oa ' (l * u> dXjdXj 


a 2 w 


Pi - £ ] 


+ e 


ikl 3 xk 


_a_ r ^1 ° "1 _ _ 1 _ . 

ax., i p o Dt 4 ax^axu F l J 


a 2 w. 


3Xj3Xj 


(4.39) 


For generality, the bracketed terms must vanish independently. Thus. 


o ur 


= 0 


(4.40a) 


D o w l ,2w l 
- - ax j ax j 


"o Dt 


- F l = 


(4.40b) 


Notice that equation (4.40b) is completely independent of w and p A , and, 
consequently can be solved separately. In fact, this is the vortical 
component of the flow, which behaves in an identical manner for both 
compressible and incompressible flows. The fundamental solution of (4.40b) 
in the non-convective form was originally developed by Ladyzhenskaya 
(196 9). This provides the basis for the development of the convective 
solution to (4.40b), as will be seen subsequently. Hcwever, next attention 
turns to the dilatational component of the flow. 

Hie velocity appearing in the linearized continuity equation (4.37a) 
can also be decomposed. As a result, 
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(4.41) 


°CPA , . S 2 W 
Dt y O 0x.3x i 


* 


since the divergence of the vortical component is zero. Combining 
appropriate derivatives of (4.40a) and (4.41), to eliminate the variable w, 
yields the following third order differential equation for pressure: 


a2 Pj , ,(U2u) D o 
3x i 8x l P 0 C 2 “= 8x i ax i c 2 Dt 2 

where 


n = o , 


n = ( 


r(x+2n) 


p o c 0 


a 2 !, 

3x^3x^ 


r o 
“ 2 Dt 


) + 


a 2 f 

3x i 3x i 


(4.42) 


(4.43) 


with c Q representing the speed of sound in the perfect gas at the reference 
state and 


r:^> l 


(4.44) 


.2 ^O 


(4.45) 


The fundamental solution of (4.42), even in the non-convective form, does 
not appear to exist in the literature, although an attempt was made 
recently to obtain the nonconvective form by Piva and Morino (1987). 

Actually, the solutions of (4.42) that are required for the boundary 
element formulation are those due to instantaneous point mass sources and 
point forces. Furthermore, in addition to the pressure response, the 
velocity field corresponding to these sources and force must be determined. 
In all cases, the results can be determined directly from the solution of 
the equation 
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*\ + Y(X+2 ») % . 

2 Dt 


3 2 p 


U 


a V x i 


p o c o 


3X.3X. 


cl Dt 2 


- 6(x-$)Mt -z) = 0 


(4.46) 


where the scalar variable P 0 is introduced along with the usual generalized 
function 6. The subscript. U, is merely a reminder that p y is a uniformly 
moving medium solution- Equation (4.46) is a scalar damped wave equation, 
which has an approximate fundamental solution of the form 


’u -;%[<«♦ ^ «c 0 t-V ] 


(4.47) 


where 




(4.48a) 


t' = t-T 

(4.48b) 

*U = (y i -U i t')(y i -U i t') 

(4.48c) 

y i - x r*i 

(4.48d) 

i = (c 2 t' 2 -R^) 1/2 . 

(4.48e) 


The presence of the Heaviside and delta functions in (4.47) establishes the 
hyperbolic nature of the dilatational response. Thus, p^ portrays the 
propagation of a scalar wave in a moving medium. Furthermore, the 
appearance of the convective radial distance in the arguments of H and 5 
leads directly to shock phenomena. As a result, equation (4.47) is 
appropriate for supersonic, as well as, subsonic flew. 

Consider, initially, the medium subjected to a unit pulse body force. 
In two-dimensions, let 

f* = 6(x-V^(t-x)e i (4.49a) 


45 


= 0 . 


(4.49b) 


From Gel'fand and and Shilov (1964), equation (4.49a) can be 
alternatively as 


? = *<t-r) a 2 

r i 2n SXjdXj 


(In r)e^. 


which, in light of (4.38b), yields 


6(t-t) 

2n 




(In r)ej 



S(t-T) 

2n 


a^ r)e j 


Then, the pressure field can be determined by using (4.49d) in (4 
subsequently, (4.42). From the result, 

\ + J"” f 6(x-?)6(t-t) ] e 

c 2 Dt 2 dx i 

and (4.46), it is evident that 




y(X+2u) _o 


Vo 


» 2 Pa 


Dt 3x i 3x i 


P A = " a^i * 

Additionally, eliminating the Laplacian operator in (4.40a), by 
(4.41), produces 


or 


V Va 

P 0 Dt * P c Dt 


Pa - £ - 0 


D w 
o w 


Dt 


r(X+2n) D qPa _ £a + 6(t-T) 
.22 Dt p- 2np_ 

p 0 c 0 F 0 ^O 


3x7 (ln r)e j 
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written 

(4.49c) 

(4.49d) 

(4.49e) 
43) and, 

= 0 , 

(4.50) 

(4.51) 
employing 

(4.52) 

(4.53) 



The solution, w. of equation (4.53) can be found by integrating over time 
within a uniformly moving media. Finally, the dilatational component of 
the velocity is determined via (4.3 8a) as 


n (dil) _ aw 
u i " ax. 


which from (4.53) can be written 


(4.54) 


(dil) = £ r<x+2|Q 


a 2 P, 


u 


2 2 
P c 
K o o 


ax i 3x j 


L ax, ax, (p u _a u )dt ] e i 


-o ax.ax.j 


(4.55) 


with 


S(t-T) 
°U ~ " 2n 


In r. 


(4.56) 


and Py defined in Appendix C. Again, the subscripts, U, signify that the 
solutions p and a should be expressed in convective coordinates. 

To complete the unit force solution, the vortical component of the 
velocity must be added to (4.55). In this case, the equation of interest 
is (4.40b) with specified by (4.49e). Thus, 


5 0 Dt 




a 2 w x 


+ e 


lij 


s(t-T) _a_ 
2n 3x, 


(In 


r)e. 


(4.57) 


Hie solution to (4.57) can be determined in terms of a scalar which is 
the fundamental solution of the convective heat equation 


D o <a u 

P 0 Dt +|1 3 Xj 3 Xj 


8(x-$)6(t-T) 


0 , 


(4.58) 


detailed in Appendix C. In particular, from (4.57) and (4.58), 


aV, 

ax i Px i 


3Uy 

"lij ax t '3 


(4.59) 
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Substituting (4.59) into (4.57), and then taking the curl of the result, 
produces 


-2 , ^ 

p o Dt le mkl 8x k 


,2 

8 My 


r 1 ~U _ a . 6 (t— c) . 1 

e mkl e lij L M ~ 3x k 8x i 1 4nr 1 J e j * 


(4.60) 


which from (4.38a) leads to the following form of the vortical component of 
the velocity 


u (vort) = f e mkl e lij f fc 8 2 
°m 




(4.61) 


Again, the subscript U is a reminder that the time integration should be 
performed for a uniformly moving media. 

To summarize, the unit instantaneous point force solution can be 
written, from (4.51), (4.55), and (4.61), in the following form 


u. 


( 


r(x+2ii) *X 


P 2 c 2 

F oo 


8x^8Xj 


rO 


3x i 3x j 


*Ptr°u* 


8 Z n 

~ e ikl e lmj ax k 8x m (tl “u" a u ) J dx )e j 


(4.62a) 



(4.62b) 

(4.62c) 


This completely defines the fundamental solutions pertaining to point 
forces, however, instantaneous point mass source solutions are also 
required. Returning to (4.43) and letting 


= 6(x-$) 6(t-T) 


(4.63a) 



leads to the simple result that 


r(X+2n) l P U _ x_ Vo 


v o 9Xl8Xi 


,2 Dt 


(4.64) 


The corresponding velocity can be determined most easily by returning to 
equations (4.40a) and (4.41), and eliminating the pressure. This produces 


3 2 w + r(X+2n) ^o , 9 2 w , 

3x i 3x i p Q c 2 ** dx i dx i 


D 2 w 

6(x-?)6(t-e) 

_2 2 2 

C 0 C 0 


0 , (4.65) 


which when compared with (4.46), establishes 


w - 1- ft 
w 2 P U 
c o 


(4.66) 


Additionally, since (4.40b) is independent of p A , w, and ^ , the vortical 
component of velocity 


3W. 

e ijk 9x7 “ 0 ’ 

and the velocity field becomes becomes 
1 c o 9Xl 


(4.67) 


( 4 . 68 ) 


Equations (4.64) and (4.68), along with 


© A - 0 , (4.69) 

comprise the instantaneous unit mass source fundamental solution. 

The final item that is required involves the response to an 
instantaneous unit heat source. In this case, 

~i) = 0 (4.70a) 
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(4.70b) 


f i ~ o 


5 = 5(X-5)6(t-r) . 


(4.70c) 


Then, from (4.37), 


P A = 0 (4.7la) 

= 0 , (4.71b) 

and © A is simply the solution to the convective heat equation. 

It is convenient, at this point to collect the fundamental point 
force, mass source, and heat source solutions into a tensor g^, where for 
the Dirac delta functions in the infinite space. 


U a ' sVk 



(4.72) 

U « ‘ ,U 1 u 2 

P 

e) T 

(4.73a) 

= { ^1 ^2 

* 

il T . 

(4.73b) 


The superscript U denotes that is a moving medium solution. 
Furthermore, the subscripts a and p vary from one to four, while in the 
following i and j vary from one to two. Additionally, the subscript p 
always takes the value three and the subscript © is four. Then, 



r o u 

g ij 

4 

_U 1 
g i© 

«S» - 

gl 

gU PP 

„U 

g pe 


. g ©j 

g ©p 

u 

9©e J 


The individual components of g^ are detailed in Appendix C. It should be 
emphasized that these are moving force and source fundamental solutions 
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and, as such, are quite involved. The explicit form of these kernels have 
recently been obtained, however the accurate numerical evaluation of the 
functions involved at high reference velocities (U^) still requires some 
additional effort. 

It may be recalled that in previous work (e.g., Dargush et al, 1987; 
Dargush and Banerjee, 1988c, 1989), all of the convective terms were brought 
to the right-hand side and included as body forces and sources. The 
corresponding fundamental solutions then involve instantaneous stationary 
point forces and sources acting in a stationary medium. These solutions 
remain time -dependent, but take a much simpler form than the convective 
Green’s functions presented in Appendix C. Unfortunately, except in the 
low to medium Reynolds number range, the stationary fundamental solutions 
do not contain enough of the physics of the problem to produce numerical 
solutions. (This will be evident in a number of examples in Section 4.7.) 
On the other hand, the convective fundamental solutions do capture the 
nature of high velocity flews, although this is not at all obvious due to 
the complicated form of the convective kernels. However, the simplified 
fundamental solution highlighted in Section 4.3.4 for convective heat 
transfer will provide sane additional insight. 

4.3.2 Incompressible Thermoviscous Flow 

In the incompressible case, the pressure becomes superfluous and is no 
longer needed as a primary variable. Additionally, the dilational 

component of the velocity vanishes. As a result, the convective 

fundamental solution for incompressible thermoviscous flow can be written 

u « - <4 - ,5 > 

where 

u o = tu x U 2 e) T (4.76a) 
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fj, . (f, l 2 J) T . 


(4.76b) 


with a and 0 varying front one to three, and the subscript & set to three. 
The kernel g^ is detailed in Appendix D. Once again, the development of 
techniques for the numerical evaluation of these kernels is still underway. 
Meanwhile, the stationary medium fundamental solutions, pertaining to 
continuous point forces and sources, are defined in Appendix E. 

4.3.3 Incompressible Viscous Flow 

Under isothermal conditions, the temperature is not required as an 
independent variable and the corresponding degree of freedom can be 
eliminated. The convective incompressible viscous flow fundamental 
solution is then equivalent to g^j from Appendix D. 

4.3.4 Convective Heat Transfer 

The final case of convective heat transfer will be examined in some 
detail. As will be seen, the fundamental solutions are manageable, yet 
still reflect several aspects of compressible thermoviscous flow. To 
begin, the reference velocity Uj^ is introduced to (4.33) to modify the 
convective derivatives. Thus, (4.33) becomes 


pc, 


V 

V Dt 


- k 


a 2 e 

9x i 3x i 


J 


where, again. 


(4.75) 


Dt 


a 
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TZ + u. 


a 

ax, 


(4.76) 


The fundamental solution, g u , due to an instantaneous point source, 
obtained from 
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D 0 g u ,20 

ut ' k ax^Xj. " 8(x_ ? )8(t " T) » 

(4.77) 

is a well-known result. A slight generalization of 

the solution presented 

in Carslaw and Jaeger (1947) produces, in two-dimensions. 

e -R u /4ct' 

g (x i ,t, S^t) - 4nk fc , 

(4.78) 

where 



(4.79a) 

t' - t-t 

(4.79b) 

Rj = (y i -U i t')(y i -U i t') 

(4.79c) 

yi - *i-h • 

(4.79d) 

The steady- state response can be obtained from 

(4.78) by integrating 

over r. Thus, 



(4.80) 

which simplifies to 


e u in /2c 

G " <x i-«i > ■ 2nk V™*' ■ 

(4.81) 

where 


R = (y iyi )1/2 

(4.82a) 

U - (U i U i ) 1/2 

(4.82b) 


and K q is the modified Bessel function of the second kind of order zero. 
It is of interest to compare (4.81) with its stationary counterpart. Of 
course, for a heat source in a stationary medium, the fundamental solution 
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is just the potential flow kernel 




(4.83) 


Figure 4.1 provides a comparison of the two kernels. G° and G, for a source 
point at the origin, lhe kernel values are plotted for field points along 
the Xj-axis, and in the convective case for a medium moving uniformly in 
the Xj-direction with a velocity of ten. Notice, in particular, that the 
static response is syirmetric about the source point, however the convective 
response is nagnified ahead of the source point, but greatly reduced behind 
it. This latter phenomenon is just the Doppler effect applied to moving 
heat sources. Thus, as illustrated for points on the positive Xj-axis in 
Figure 4.1, the strength of an oncoming source appears to be intensified. 
On the other hand, the source has already passed the points on the negative 
x^axis, and a quick silencing is apparent. 

Interestingly, from another vantage point, the convective Green's 
function G u can be viewed as the boundary element counterpart of the so- 
called 'upwinding' techniques that are required in finite difference and 
finite element approaches to convective problems. The distinguishing 
feature is that G u embodies an analytical form of upwinding, while the 
other two methods use ad hoc representations. As a result, a boundary 
element formulation based upon G u will have a significant advantage for 
convection-dominated problems. 

The transient convective diffusion kernel can also be formed by 
integrating (4.78), but this time from zero to t. The result is a two- 
dimensional fundamental solution, which can be written in series form as. 


G u (x i ,t; ^, 0 ) 


e U iyi /2c - { - VjPi t/4c) n 

4nk * rfi E n+l 

n=0 


R 
4ct 


(4.84) 
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where E n+1 is the exponential integral of order n+1. Figure 4.2 compares 
the steady-state kernel (4.81) with this transient kernel for several 
values of t. Note that the Doppler effect is still quite pronounced. 

Before closing this section on fundamental solutions, it should be 
emphasized that behavior similar to that displayed in (4.81) and (4.84) 
will be included in the convective thermo viscous kernels, since the scalar 
Green's function G u provides the basis for the development of the more 
complicated fundamental solutions. In fact, for the incompressible 
theories, G u is the only scalar Green's function that is needed. (More 
precisely, a change in material constants is required to produce of 
equation (4.58) from G U .) However, for compressible flow a second scalar 
fundamental solution, comes into play for the dilatational component of 
the flow. As mentioned previously, this latter solution involves the 
propagation of a damped wave, which at high velocities produces shock 
phenomena. 


4.4 Integral Representations 

4.4.1 Compressible Thermoviscous Flow 

The desired integral representation for general compressible 
thermoviscous flow can be derived directly from the set of governing 
differential equations. First, however, a convenient differential operator 
notation is introduced. As a result, equations (4.37) are rewritten as 


L aP u p + 


f a = 


where, again 


Up = (Uj u 2 p $} T 

f„ - {?! ? 2 ; Ji T 


(4.85) 


(4.73a) 

(4.73b) 
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(4.87C) 
(4. 87d) 


&--K 


(4.87e) 




(4. 87f ) 


T U = 0 

L ej 0 


(4. 87g) 


i” « 0 
0p 


(4.87h) 


rU 


ee " _p o c v Dt + k ax^Xjp 


(4. 87 i ) 


Then, using L^p to operate on the fundamental solution g^ ft of (4.74) 


1 ap 


produces 


L “a9p r + - o . 


(4.88) 
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In (4.88), the subscript y also varies from one to four and Kronecker's 
delta function has been generalized in an obvious manner. 

(the governing equations (4.84) must, of course, hold for all points of 
the flow region at every instant of time. Therefore, the lefthand side of 
(4.85) multiplied by an arbitrary function g a<y , and integrated over time 
and space must remain equal to zero. That is, 

t 

< V l “o u k * h > - J 0 J v 5«r<^v I . ,d ' , ‘ h • 0 • <4 - 89) 

where the standard notation for the inner product of two functions has been 
introduced. Returning to the explicit forms of the differential operators, 
this becomes 

l 0 J v '9ir '-“o' 1 ! - OoVi.m + a+ i‘ )u n.ta * “ u i,mn - “p,l * ? 1> 

+ ^pr ^"Po^.m ^ ” ^V^.m + ^p* 

+ [ -Po c v u o - p o c v u m u 0,m + ku e.im, + J e ]dvdT “ 0 * (4 * 90) 

in which commas represent spatial derivatives and superposed dots are 
partial derivatives with respect to time. Next, the divergence theorem can 
be applied, repeatedly, to the applicable terms in (4.90) to transfer 
spatial, as well as, temporal derivatives from Up to g Qiy . As a result, 
equations (4.90) are transformed into 

l 0 i s '5„ r t„ - ?„ T “« ,asdT * J 0 f v <S„ T J „!<w<fc - J v I 5oY u £ , o lav 

+ U/'^Y + p o°Tn5iY,m + * ^ir.rrm + PoV,.i lu i 

+ ^my,m + ^py + ^m^py.m^ 
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(4.91) 


+ ^ p o c v^©r + p o c v^m^0Y.m + = 0 ' 

where 


fc i ~ Xu m,m n i + p ^ u m,i +u i,m )n m " 

(4.92a) 
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4 

D 
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It 

(4.92b) 

fc © = ku «.nftn " VvW 

(4.92c) 

^iY - ^W.m n i + 2p ^iY,m n m + Po^pY n i 

(4.93a) 

? PY = 5 iy n i 

(4.93b) 

^©y = ^©Y.nftn 

(4.93c) 

U £ = {p O U l p O U 2 °p p O C V U © )T 

(4.94) 


with defined as the unit normal to the surface S at X. To complete the 
derivation of the integral equation for any point $ interior to S at time 
t, the last volume integral appearing in (4.91) must be reduced to 
“ u y ($»t). Hiis is accomplished, if 

< ^aY'^ > “ -V5. fc > (4 * 95) 
or after making use of the properties of the delta function 


*?«9ar + 8 p/ (x ‘ 5)8(t ' T) " 0 * 

where the operator has components 

■ *ji»o K * <w> * ‘ji 11 5^ 

*?p - p o Ttq 


(4.96) 


(4.97a) 


(4.97b) 
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= 0 

(4.97C) 
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= 0x i 

(4. 97d) 
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L pe 

= 0 

(4.97f) 
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= 0 

(4,97g) 
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= 0 

(4.97h) 
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L ee 

V a 2 

0,3 

■ P 0 C V Dt + K * 

(4. 97i) 


Formally, £^ a is called the adjoint of the original compressible 
thermoviscous differential operator and g ay defined by (4.95) is the 

adjoint Green's function. This adjoint Green's function can be obtained 
simply by suitably transposing the fundamental solution presented in 
Section 4.3.1. That is, 

g ay (X,t;5.r) = g ya (5.r;X,t) . (4.98) 

Substituting (4.98) into (4.91) produces the desired integral equation, 

u r u.t> - | s (g“ a *t 0 -f“ a *u a )ds * f v !g“ a -E 0 iav (4.99) 

in which, for simplicity, the initial conditions have been assumed zero. 
The • in (4.99) once again symbolizes a Riemann convolution integral. 

Notice that this integral equation for compressible thermoviscous flow 
has a similar form to that for thermoelasticity as shown in equation (3.2). 
However, in (4.99), a volume integral is retained to include, in 
particular, the nonlinear body force terms. 
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A derivation of the integral representation for the incompressible 

flow theories would follow the same lines as that just presented, and 

therefore, will not be repeated. In fact, a generalized integral equation 

identical to (4.99) would result. The only differences are in the explicit 

form of the fundamental solutions g^ a and in the corresponding definitions 

of the functions f u and t . 

yet o 

As may be recalled from (4.3 5), a portion of the convective effects 
are included in the body forces f a . Assuming for the moment that this is 
the only non-zero component of f Q , then the volume integral in (4.99) can 
be rewritten as 




(4.100) 


Applying the divergence theorem to the right-hand side of (4.100) produces 

jyfeVwjVj]*' ' l s [ g ?«' CU a U j n jl dS - J v [9?.,j*<’ u « u j] av • (4 - 101) 

since, for the incompressible case, u..^ is identically zero. Finally, 
equation (4.99) becomes 




(4.102) 


where 


t ' = t - pu u-n- . 
a a r a ] j 


(4.103) 


Notice, in particular, that (4.102) no longer involves velocity gradients. 
Consequently, from a computational standpoint, (4.102) is an attractive 
alternative to (4.99). 


A similar integral formulation can also be developed by utilizing the 

stationary medium fundamental solutions g . in this case, the reference 

ya 
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velocity IK may still be used, but now the entire convective derivative 
must be included in the body forces. As a result, the integral equation is 
written 


Vi.t) = i s [g 7a -t;'-f Ta -u o ]as * J v [9 ra . j *pu <1 v j ]» 
in which 

fc i* ■ l a - 0 U a V j n j 
and Vj is the total velocity. 

4.4.3 C o nvec t i ve Heat , T r a nsfer 

In this simplest case, equations (4.99) reduce to 

0(?,t) * f [-g%+f U *e]dS + f [g U *J]dV (4.106) 

s J v 

where g u is defined by (4.78) and 



(4.107) 


(4.108) 


Meanwhile, under steady conditions, equation (4.106) simplifies to 

©U.t) = f [-G U q+F U 0JdS + f (G U 3JdV (4.109) 

5 J V 

in which 


cU 


with G° 


given by (4.81). 


(4.110) 


4.5 Numerical Inplementation 

The numerical treatment of the equations in thermoviscous fluid 
dynamics follows very closely that described in Section 3 for transient 
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thermal stress analysis. However, now due to the volume integral appearing 
in (4.99), (4.102) or (4.104), the interior must be subdivided into cells. 
The geometry of each cell is defined by nodal points and quadratic shape 
functions. In two-dimensions, six and eight-noded cells are available. 
Meanwhile, either a linear or quadratic variation can be employed for the 
functional representation. Details of the techniques used for cell 
integration can be found in Mustoe (1984). 

Just as for the thermoelastic case, a set of algebraic equations can 
be developed by writing the integral equation at each global node. 
However, now interior, as well as, boundary nodes must be included, and the 
resulting equations become highly nonlinear due to the convective terms. 
After the collocation process is complete, the final system of equations 
can be expressed in matrix form as 

g b (x,u) = A b x + I^o 0 - <Pt° - B b y = 0 (4.iila) 

for boundary points, and as 

g U (x,u) = U + A U x + D U a° - G U t° - B'V = 0 (4.111b) 

for interior cell points, where the vectors <x° and t° have components 
defined by 

"ip ‘ » u i u o 

i-O O „ 

" a ip n i 

at each boundary and interior point. Once again x and y are the known and 
unknown boundary quantities, while u is the interior velocity vector, and 
the matrices A, B, D and G are developed from the integrals of the kernel 
functions appearing in (4.99), (4.102) or (4.104). At present, only 
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(4.104) has been implemented as a segment of the general purpose boundary 
element program, GP-BEST. 

Initially, an iterative algorithm, along the lines of those used for 
BEM elastoplasticity, was enployed to solve (4.111). However, convergence 
is usually achieved only at low Reynolds number. More generally, when 
employing the stationary fundamental solutions, the interior equations must 
be brought into the system matrix along with the boundary equations, and a 
full or modified Newton-Raphson algorithm must be utilized to obtain 
solutions at moderate or high Reynolds number. Symbolically, at each 
iteration m, 

(Ax) m 
(Au) m 

where 

3 ^ n+1 = X 111 + (Ax) 1 " 

u m+1 . U 111 + (Au) m 

and the derivatives on the lefthand side of (4.112) are evaluated at 
(xV). In the numerical implementation, the above equations are arranged 
to form a block banded system matrix for efficient multi-region solutions. 

It is anticipated that once the convective viscous kernels are 
implemented somewhat different solution strategies will be more 
appropriate. For example, at high velocities the system matrix will become 
sparse. In that case, bandwidth minimization is required and iterative 
equation solvers become quite attractive. 
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4.6 Coupling of Solid and Fluid 

The coupling of the solid and fluid phases is most readily 
accommodated via the concept of the generic modeling region. Thus, the 
fluid-structure interface is nothing more than a boundary between two 
GMR's. In the simplest case, temperature, flux, and tractions are matched 
across the fluid-structure interface, while a temporal approximation is 
introduced to relate boundary displacements of the solid to the 
corresponding fluid velocities. However, additional sophistication is 
possible. For example, thermal resistance can be introduced to model the 
effects of coatings. 

4.7 Examples 

4.7.1 Parallel Flw 

The two-dimensional parallel flow in a duct is a good verification 
problem for incompressible computational fluid dynamics codes. It has a 
simple analytical solution which can be used to test many aspects of 
programs. The convective terms disappear in the nonlinear solution, hence 
linear and nonlinear velocity profiles should be identical (Tadmor and 
Gogos, 1977). 

As an example of a typical version of this problem. Figure 4.3 
illustrates a 10 cell mesh with two regions. This simulates a plate 
sliding along the top of the fluid in pure shear. Pure shear tractions are 
applied at inlet and exit. Viscosity is unity and density is incremented 
to increase the effect of the convective terms in the equations. Newton 
iteration is used to converge to the nonlinear solution. It should be 
noted that this problem does not require this degree of refinement. This 
model merely tests many aspects of the computer program. 

Figure 4.4 illustrates the linear velocity profile at the exit of the 
region. For density below iooo the linear profile is reproduced exactly. 
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4.7.2 Priven Cavity 

The two-dimensional driven cavity has become the standard test problem 
for incompressible computational fluid dynamics codes. In a way, this is 
unfortunate because of the ambiguities in the specification of the boundary 
conditions. However, numerous results are available for comparison 
purposes. 

The incompressible fluid of uniform viscosity is confined within a 
unit square region. The fluid velocities on the left, right and bottom 
sides are fixed at zero, while a uniform non-zero velocity is specified in 
the x-direction along the top edge. Thus, in the top corners, the x- 
velocity is not clearly defined. To alleviate this difficulty in the 
present analysis, the magnitude of this velocity component is tapered to 
zero at the corners. 

Results are presented for the 144 cell boundary element model shown in 
Figure 4.5. Notice that a higher level of refinement is used near the 
edges. Spatial plots of the resulting velocity vectors are displayed in 
Figures 4.6, 4.7, and 4.8 for Reynolds numbers (Re) of 100, 400 and 1000, 
respectively. Notice that, in particular, the shift of the vortical center 
follows that described by Burggraf (1966) in his classic paper. A more 
quantitative examination of the results can be found in Figure 4.9, where 
the horizontal velocities on the vertical centerline obtained from the 
present analysis (i.e., GP-BEST results) are compared to those of Ghia et 
al. (1982). It is assumed that the latter solutions are quite accurate 
since the authors employed a 129 by 12 9 finite difference grid. It is 
apparent, from the figure, that the present boundary element model has some 
difficulty in capturing the sharp knee of the curve at Re = 400. This 
becomes accentuated as the Reynolds number increases, and consequently, a 
finer mesh is required. It should be noted that the simple iterative 
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algorithm fails to converge much beyond Re = 100. Beyond that range the 
use of a Newton- Raphson type algorithm is imperative. 

In order to obtain more accurate solutions at higher Reynolds number, 
the refined four region 324-cell boundary element model shown in Figure 
4.10 was also analyzed. This provides a significant improvement in the 
results. For example, at Re = 1000, as seen from Figure 4.11, the 
secondary vortex in the lower right-hand corner is clearly visible. 
Additionally, the resulting horizontal velocities are compared to Ghia et 
al (1982) in Figure 4.12. Now, even the solution at Re=l000 is in 
excellent agreement. 

4.7.3 Converging Channel 

The two-dimensional incompressible flew through a converging channel 
also possesses a well known analytical solution which is purely radial 
(Millsaps and Pohl hausen, 1953). A comprehensive finite element study of 
this problem has been made by Gartling, et al (1977). 

The boundary element model is shown in Figure 4.13. The mesh contains 
96 cells and is divided into two regions. The boundary conditions were 
modeled using an exact specification of the boundary conditions appearing 
in the analytical solution (Fig. 4.13). Viscosity is unity and tractions 
and density are incremented to reach higher Reynolds numbers. The Reynolds 
number for this problem is defined as 

pRjV,(R.) 

R - — ^ — i - 

e n 

where Vj(r^) is the maximum velocity in the region, which is -24.0 for the 
problem solved here. 

Figure 4.14 illustrates the results for two Reynolds numbers, 
indicating good accuracy along the entire width of the channel. Not only 
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are the velocities accurate, but the pressures and tractions are very 
accurate also. 

It has been observed that finite element versions of this problem have 
several peculiarities which prevent the analytical solution from being 
reproduced. First of all, velocities are often specified at the inlet and 
at the wall and centerline, ambiguous boundary condition specification 
results. Also, typically a parabolic "fully developed" velocity profile 
is often specified at the inlet. However, the nonlinear solution has a 
flattened velocity distribution across the width of the channel (see Fig. 
4.14). Hence, the analytical solution cannot be reproduced exactly if the 
"fully developed" profile is specified at the inlet. Also, the finite 
element modelers of this problem usually leave out the traction 
distribution at the exit and specify zero tractions there. This also gives 
rise to non-radial flow. 

The reason for so much interest in the converging flow problem is that 
it is one of the few problems possessing an analytical solution. However, 
by specifying a model which does not correspond to this problem, as in the 
finite element case, one cannot accurately compare results to the 
analytical solution. Any such comparisons are merely qualitative. In this 
light, the boundary element model here has utilized an exact model of the 
boundary conditions appearing in the analytical solution. This way an 
accurate and meaningful comparison can be made. 

4.7.4 Flow Over a Cylinder 

Next, an example of unconfined flow around an obstacle is considered. 
In particular, the oft-studied case of a unit diameter circular cylinder is 
examined. The boundary element mesh is illustrated in Figure 4.15. Notice 
that three distinct regions are evident. The smallest region, labelled 
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GMRl, represents a thennoelastic thick-walled cylinder. Only the surface 
of the solid is discretized. The next region, GMR2, models a thermo viscous 
fluid in the vicinity of the cylinder. In GMR2 voltane cells are required 
due to convective body forces. However, sufficiently remote from the 
cylinder, these body forces become negligible and once again a boundary- 
only region, in this case GMR3, is valid. 

Steady-state velocity vector plots are displayed in Figures 4.16 and 
4.17 for Re * 20 and 40, respectively. The recirculating zone, behind the 
cylinder, is clearly visible. 

Additionally, the problem was extended to include thermal effects. 
The temperature of the fluid at inlet was specified as looo°C, while that 
at the inner surface of the hollow cylinder was maintained at 0°c. The 
effective heat transfer coefficient between the fluid and solid can then be 
obtained from the resulting temperature and flux at the outer surface of 
the cylinder. The distribution of the nondimensional Nusselt number (Nu) 
around the circumference is plotted in Figure 4.18. These curves agree, at 
least, qualitatively with the experimental results of Eckert and Soehngen 
(1952). Of course, if the purpose of the analysis is to determine the 
temperature and stress in the solid, then there is really no need to 
compute the heat transfer coefficients. The desired solid temperatures and 
stresses come directly out of the analysis. 

4.7.5 Flow Over an Airfoil 

As a final example, the themoviscous flew over a NACA 0018 airfoil is 
considered. The boundary element model shown in Figure 4.19 once again 
utilizes symmetry and employs the multiregion concept with cells confined 
to the vicinity of the airfoil. The airfoil is heated externally by a hot 
gas, flowing from left to right, at unit temperature, and cooled to zero cm 
the surface of an internal cooling hole. The conductivity of the airfoil 
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is set to one hundred times that of the hot gas, while unit values are 
assumed for the fluid density and viscosity. 

The resulting steady-state velocity distribution at Re = 150 is 
displayed in Figure 4.20, while Figure 4.21 details the velocity profile 
just ahead of the leading edge of the blade. It should be noted that lift 
and drag can easily be calculated, since during the analysis, tractions are 
determined all along the blade surface. Next, shaded temperature contour 
plots of the region surrounding the airfoil are presented for Re = 10 and 
150 in Figure 4.22. In the latter diagram, the hot regions are black, 
while lower temperature locations appear white. The effects of convection 
are visible downstream of the airfoil. Lastly, the surface temperature of 
the airfoil is plotted in Figure 4.23. Notice that the overall temperature 
increases with Reynolds number. In this particular case, the distribution 
is strongly influenced by the location of the single internal cooling hole. 

When the Reynolds number is elevated further, the convective terms 
begin to dominate. In this flow regime, the physics of the problem demands 
that convective effects must be incorporated in the kernel functions. This 
is, in fact, true for all of the viscous flow examples presented thus far. 
As mentioned earlier, the inclusion of convection in the kernel functions 
is analogous to the upwinding techniques that are required in finite 
difference and finite element analyses. 

The development and numerical verification of these convective 
thermo viscous flow kernels is new underway. However, the thermal portion 
of the new kernels, detailed in Section 4.3.4 and 4.4.3, has been 
implemented and provides some interesting results. 

As an illustrative example, a convective heat transfer analysis was 
conducted for a pair of NACA 0018 airfoils in a uniform flow field. The 
boundary element model of the airfoils is shown in Figure 4.24. The hot 
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fluid once again flews from left to right, while the airfoils are cooled on 
their inner surfaces. It should be emphasized that with the assumption of 
a uniform fluid velocity, the problem permits a boundary-only analysis. 
Thus, the only mesh that is needed is that displayed, in Figure 4.24. 
However, a number of interior points were added in the flow field for post- 
processing purposes. 

Figure 4.25 depicts the temperature distribution in the fluid 
surrounding the airfoils at a Feclet (Pe) number of ten, where 



with fluid velocity U, chord length 1, and thermal diffusivity of the fluid 
c. Meanwhile, Figures 4.26 and 4.27 present the temperature field for 
Pe=lOO and 1000, respectively. Strong convective effects are evident at 
the higher Peclet numbers. Finally, in Figures 4.28 and 4.29 the angle of 
attack is modified to 10° and 20° while maintaining Pe=l000. 

It should be reiterated that the results shown in Figures 4.25-4.29 
are based on a uniform flow field. Thus, the effects of viscosity have 
been ignored. However, the new convective thermo viscous kernels, when they 
are available, will have the same character as those for convective heat 
transfer, and hence, should provide a means for obtaining accurate high 
velocity solutions. 
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FIGURE 4.5 


DRIVEN CRVITY 
Boundary Element Model 



FIGURE 4.6 
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Steady-State Re “ 100 
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FIGURE 4.7 







HORIZONTRL VELOCITY at X=0.5 







FIGURE 4.10 


DRIVEN CAVITY - FOUR REGION MODEL 
Boundary Element Model 



DRIVEN CAVITY - FOUR REGION MODEL 
Lower Right Corner at Re = 1000 
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HORIZONTRL VELOCITY at X=0.5 





FIGURE 4.13 


CONVERGING CHANNEL - PROBLEM DEFINITION 
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FIGURE 4.16 
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FIGURE 4.17 
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FIGURE 4.25 
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CONVECTIVE THERMAL FLOW ( PE =10, ANGLE 



FLOW PAST A PAIR OF AIRFOILS 
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CONVECTIVE THERMAL FLOW ( PE =4000, ANGLE 
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CONVECTIVE THERMAL FLOW < PE =1000, ANGLE 


5. STOMUff 


Significant advancements have been made in the last twelve months 
toward the development of an integrated boundary element method for hot 
fluid-structure interaction. For the solids portion of the problem, the 
formulation is well developed. The boundary-only time domain thermoelastic 
formulation, detailed in Section 3 and Appendix B, was completed in the 
previous year. However, a number of enhancements have been incorporated to 
make the numerical implementation more efficient, more accurate, and to 
increase its applicability. For example, regarding computational aspects, 
full advantage is now taken of the uncoupled nature of the thermoelastic 
theory, so that convolution is only carried out on the temperature and flux 
related quantities. Additionally, for time steps beyond the first, a much 
reduced level of numerical integration is employed to evaluate the 
completely non-singular kernel functions. Meanwhile, extensions of the 
basic formulations have been made to include several practical facilities, 
such as time-dependent ambient temperatures, thermal resistance between 
regions to simulate coatings and air gaps, and the introduction of region- 
by-region reference temperatures. The resulting code has also gone through 
another round of verification testing, which has greatly improved its 
reliability. 

The primary emphasis of the work performed under this grant has, of 
course, been directed toward the fluid, since boundary element applications 
to fluids are at a much less developed state. Considerable progress has 
been made on two fronts. Hie first major area involves improvements and 
extensions of the incompressible thermoviscous formulation originally 
developed last year. During the past twelve months, the accuracy and 
efficiency of the numerical integration has been significantly upgraded, 
the volume integrals have been rewritten to eliminate the need for 
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computation of velocity gradients, reference velocities and temperatures 
have been introduced on a region-by-region basis, and a Newton-Raphson 
algorithm has been developed to solve the highly nonlinear set of 
equations. . Hie result, as is evident from the examples of Section 4, is an 
accurate general purpose boundary element approach to problems of 
thermoviscous flow in the low to medium Reynolds number range. As such, 
this development represents the first of its kind for this class of 
problems. 

However, during the course of this work, it also became evident that 
the stationary media fundamental solutions of Appendix E do not contain 
enough of the physics of the problem at high Reynolds number. Moving media 
fundamental solutions and integral formulations are imperative for higher 
speed flows. Since these fundamental solutions do not exist in the 
literature, considerable effort has been expended toward their derivation. 
Approximate forms have been obtained for oompressible thermoviscous flew, 
and are presented in Section 4.3. It should be emphasized that these 
convective solutions contain an analytical representation of upwinding and, 
for compressible flow, shock. The development of techniques for the 
numerical evaluation of the convective kernels is now underway. Meanwhile, 
the thermal portion pertaining to convective heat transfer, in a knewn flow 
field, has been completely implemented, lhis new formulation not only has 
produced some interesting results, but also provides considerable optimism 
for the success of the convective media approach to high speed 
thermoviscous flow. 
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6. WORKPLAN FOR THE NEXT YEAR 


Based upon the experiences of this past year, future emphasis will 

naturally be placed upon the convective media approach, although some 

ongoing work on the transient stationary media algorithm will be completed. 

The following rather ambitious set of tasks are planned, in approximate 

chronological order, for the period November 1988 to November 1989: 

1. Complete development of numerical techniques (e.g., rational 
approximations, series representations) for the evaluation of the 
convective compressible thermoviscous kernels. 

2. Implement and validate the transient convective heat transfer 
formulation. 

3. Complete the investigation of transient incompressible flew using the 
stationary media approach. 

4. Implement and validate the new convective incompressible flow kernels. 

5. Develop more efficient solution algorithms (e.g., iteration methods) 
and integration schemes for high Re flow. 

6 . Implement and validate convective compressible flow kernels. 
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APPENDIX B - KERNELS FOR THEEMOEUSTICnY 

This appendix contains the detailed presentations of all the kernel 
functions utilized in the formulations contained in Section 3. lV/o- 
dimensional (plane strain) kernels are provided, based upon continuous 
source and force fundamental solutions. For time -dependent uncoupled 
quasistatic thermoelasticity the following relationships must be used to 
determine the proper form of the functions required in the boundary element 
discretization. That is, 

G” p (x-?) - G op (X-?,nAt) for n=l 

g£ p (X-$) * G ap (X-?,nAt) - G a p(X-$. (n-l)At) for n>l . 

with similar expressions holding for all the remaining kernels. In the 
specification of these kernels below, the arguments (X-$,t) are assumed. 
Ihe indices 

i,j,k,l vary from 1 to d 

a,p vary from 1 to (d+1) 

e equals d+l 

where d is the dimensionality of the problem. Additionally, 

coordinates of integration point 
5^ coordinates of field point 

y i = Xj.-q r 2 * y^ . 

For the displacement kernel, 
l 1 

G. . = - — [ (-V 1 ) - (6. •) (3-4v)ln r ] 

13 8 n n(l-v) r 2 13 

G w ■ 0 
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y-s - 


G ej ■ b ( Eofer ) 1 , F a> 94 ( " > 1 


G M - 2n <*> ' 95<1> 1 


whereas, for the traction kernel. 


. j_ _i_ t ( 2 yjyjyk n k ) _ 


ij 4nr (1-v) 


y i n i 

+ (-H 1 ) 


(l-2v) ] 


F iO = 0 




i_ ,_L_x r , y i y k n k . - 
4n X+2n^ r 2 ^ 


6 (n) - (nj)f 7 (T,) ] 


ee 


l y k n k - 

" 2^ [ HH f 8<*> 1 


In the above. 


tj = 


(ct) 


1/2 


c m Js_ 


E i <*> ■^ s r dx 


h x ( n ) = ^ (l-e" 11 /4 ) 


(l-2v) 
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h^n) E,( 4 ) 
g 4 (T,) = 2 + -V~ 


g 5 (t)) = — 


Ei(J) 


f 6 (r\) « h^n) 


f 7 ( T l ) 


tyi) ) 


fg(ri) = e 


-T1 2 /4 


For the interior stress kernels, 


E pij 1-2 v 8 ij a^ + 4 + p6 ij G pe 


2 . 3G J1 + ( f£g i + Si) . 


D 0ij 


^ •« *£f> ' M„ P P9 


where 


aG li . _i_ _i . ( ^££k . ij£i . fi£i, 

9? k 8nr n(l-v) 1 j.3 r r 


+ (-4^) {(3-4v) ] 


9? k 4n 'ka+2(i) 


41 - - 1 ( ^ j r ) < h i> - < 8 jk> <r + r } ] 
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8F n j i_ . 4 yiyiykyi n i y^k 

8 <k '4« J <>-»> 1 r < r * t * 


8 tky-iyi"i l7 ,5fWi ; ytyk n i . , , 

t 2 '* J 2 6 ij"k ^2 " s ik n j )f 2 , ' |) 


* <-^2 - Vi 1 ') 1 ' 1 1 


3F 


ei i 




a? k 4nr X+2(i' 1 ' r 3 


y k n i y i n k 8 ik y l n l - 

- <-p + -p + ■ 1 Y ~ 1){h i } 1 * 


fjh) = 2 

f 2 (tj) = l-2v 
f 3 (t| ) = l-2v 
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APPENDIX C - FUNDAMENTAL SOLUTION FOR CONVECTIVE COMPRESSIBLE THERMOVISCOOS 
FLOW 




2 

r(x+ 2 h) 9 p u 1 r fc r a 2 , ft . 

>j c j 8 x i ax j p o J o 1 9 v x j >u '°° 


8 2 n 

" e ikl e imj I st 


& 


_J 5 n 

8 Xj 


9 ej ' 0 


o . i_ !!c 
3l P c 2 8x t 


u _ y(X+2n) 9 _ y_ DpPu 


3 PP 2 ax^Xi 

w CTO A A 


,2 Dt 


9ep = 0 


9 ie “ 0 


9 pe “ 0 


9 ee " «u 


where 


°U 


S(t— c) 
2 n 


In r 


*^J ~ 4nn 


# -ty 4 c'f 

t 5 
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e-V 4ct ' 


©, 


0 4nk t' 


Q 

P U - 7T [ tl * ?T H(c 0 f-Ru) - \ 6<c 0 t 


2nR 


C' 

P. 


p o c v 


= X+2u 


t' = t-x 


bJ - (y i -o 1 f)(y l -O t f 


y t - 


? - c 2 c 2 - F* 


’‘V ] 
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APPENDIX D - FTOBftMENEAL SOUDTION PCJR CONVECTIVE 
7BEROVISGOOS FLOW 



e ikl e lmj 




9ej = 0 

_u 


g ee 


u 


where 


°U 

% 


S(t-t) 

2n 


In r 


j .'V*'*' 

4n(i t' 

j e'V**' 

“ 4nk t' 


c' = * 

p o 

t' = t-T 

R* = (yi-Uit'Xy. -t^f) 

y i = x i‘ 5 i 


p o c v 


INCOMPRESSIBLE 
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APPENDIX E - KEFNELS FOR STATIONARY HKXMPRESSIBIf THEHHCWISOOOS FLOW 


This appendix contains details of the time-dependent incompressible 
kernels, based upon stationary media, necessary for the integral 
formulations of Section 4. Notation is consistent with that defined in 


Appendix B. 

For the generalized velocity kernels. 


l r y i y i s l (e) 

G ij - 7^ [ - <6 ij> <— 


=>£>] 


G i« = ° 


G 0j = 0 


G ee " 2n ( k ) [ 




whereas, for the generalized traction kernel. 


p ij - 2^ [ (s r^* Z,4) * ts r e- t2,4 i 


- (- 


A (H(t) . 8i , . (JSl - e - 2 ' 4 , ] 


f « - 0 


r «j - ° 


'ee 2«r I ' r ’ e J * 


In the above, 


E = 


(c't) 


1/2 
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c' = ^ 


S.U) = ~ (l-e‘ e /4 ) 


(ct) 


pc. 


1/2 


z 


Meanwhile, for the interior strain rates. 


E ijk - ^ «•!> - < !± F 1 > •*!> * ' 2 n-«-' 2/4 > 

♦ t2e"* 2/4 -Sj) ] 


ijk 


1 r ,¥WA,: , 6 llWm , W,- , 

2„r 2 l ’ 1* 91 ‘ S r 2 >92 * r 2 >9j 


+ , Y i y i n k + S&jtuh + ilt^|W4n )gi 


( 'ij n k * 5 lk n j ) 9 2 - <6 jk n i>9 3 ] 


where 


gj = 4Sj - 2e 


-e 2 / 4 


g 2 = -Sj + e 


-e 2 /4 


109 



i cp i cn i tn 


l 

l 


g 3 « H(t) - Sj 

2 2 

j = 24Sj - I6e -T5 /4 - e 2 e“ e ^ 4 

2 * -4 Sl + re" e2/4 + — e" e2/4 

3 = 4s x - 2e -8 /4 - 2H(t) . 


110 


